5 de Mayo de 2021

Objetivo

En este taller estudiaremos información geográfica, su análisis y visualización. En particular

  • Aprenderemos como visualizar mapas, shapefiles, etc. en R
  • Utilizaremos datos geolocalizados para estimar la residencia de usuarios
  • Utilizaremos datos geolocalizados para estimar la ley de la gravedad de retail
  • Utilizaremos esos modelos para determinar las zonas de atracción de centros comerciales

Visualización de información geográfica

Hya muchas librerías de visualización de información geográfica. A estos sistemas se los conoce como GIS (Geografical Information Systems).

  • Una lista de estas plataformas/sistemas se puede encontrar en la Página en Wikipedia sobre GIS

  • Ejemplos importantes de estas plataformas GIS son:

    • Los productos Esri, que incluyen el famoso ArcGIS (desktop y online versions)
    • El servicio (español) Carto, un servicio online muy sencillo para mostrar información geográfica.

  • En R hay muchas librerías para visualizar y analizar datos geográficos.

    • Una lista de paquetes se puede encontrar en la Task View for Spatial Statistics
    • Principalmente hay dos categorías
      • Paquetes para tratar datos tipo (lat, lon) como sp, sf,maps, rgdal, etc.
      • Paquetes que ofrecen integración APIs de mapas como Google Maps o OpenStreetMaps: RgoogleMaps, ggmap, leaflet, etc.

  • Lo mismo en Python

Tipos de datos geográficos

Recordamos que hay dos tipos de datos geográficos

  • Raster o imágenes
    • Como por ejemplo mapas de carreteras, fotografías aéreas, etc.
  • Puntos, líneas o polígonos
    • Como por ejemplo las líneas que definen una carretera, el polígono que delimita un municipio o un punto que determina desde dónde se ha geo-localizado un tweet.

Vamos a ver cómo usar ambos.

Tipos de datos geográficos

Usando raster data

En R hay muchos paquetes para descargar, tratar y analizar imagenes geográficas tipo raster. Empecemos con los que se descargan imágenes.

Básicamente estos paquetes descargan una imagen (raster) de un mapa y permiten mostrar sobre estos mapas puntos, líneas, polígonos, etc.

Antes de pintar el mapa, vamos a bajarnos algunos puntos geolocalizados. En concreto los centros comerciales en la comunidad de Madrid (http://www.madrid.org/nomecalles/)

centros <- read.table("./data/centros_comerciales.csv",header=T,stringsAsFactors = F)
head(centros,2)
##   CMUN     ETIQUETA         MUNICIPIO
## 1    5 Leroy Merlin Alcalá de Henares
## 2    7 Leroy Merlin          Alcorcón
##                                                                          DIRECCIÓN
## 1        Ctra Nacional II, km. 34 - Calle Tamajón s/n (Centro Comercial La Dehesa)
## 2 Avenida de Europa, 26 (Centro Comercial Parque Oeste - Ctra. Nacional V, km. 14)
##         lon      lat
## 1 -3.327188 40.50292
## 2 -3.849481 40.34330

Como vemos los centros comerciales vienen georeferenciados con coordenadas lat lon

Paquetes con integración con APIs de mapas

Pintamos el mapa de Madrid con los puntos de los centros comerciales. Primero lo hacemos con el paquete RgoogleMaps que toma los mapas de Google

require(RgoogleMaps)
MadridMap <- GetMap(center = c(40.415364,-3.707398), zoom = 11,maptype="mobile")
tmp <- PlotOnStaticMap(MadridMap,lon=centros$lon,lat=centros$lat,pch=20,cex=2,col=factor(centros$ETIQUETA))

Paquetes con integración con APIs de mapas

Podemos cambiar el tipo de mapa para mostrar diferentes tipos de imágenes `maptype = c(“roadmap”, “mobile”, “satellite”, “terrain”, “hybrid”, “mapmaker-roadmap”, “mapmaker-hybrid”)

MadridMap <- GetMap(center = c(40.415364,-3.707398), zoom = 11,maptype="hybrid")
tmp <- PlotOnStaticMap(MadridMap,lon=centros$lon,lat=centros$lat,pch=20,cex=2,col=factor(centros$ETIQUETA))

Paquetes con integración con APIs de mapas

Lo mismo lo podemos hacer con el paquete ggmap (https://github.com/dkahle/ggmap). Y un poco más fácil.

require(ggmap)
madrid_bb <- c(-3.933449,40.290525,-3.476830,40.538069)
mm <- get_stamenmap(madrid_bb,zoom=11)
ggmap(mm)+geom_point(data=centros, aes(x=lon,y=lat,col=ETIQUETA))

Paquetes con integración con APIs de mapas

Podemos hacer esto con visualización interactiva a través del paquete leaflet https://rstudio.github.io/leaflet/ Los mapas por defecto son de OpenStreetMap

m <- leaflet(data=centros) %>% 
  addTiles() %>% 
  addMarkers(~lon, ~lat, popup = ~ as.character(ETIQUETA))
m

Paquetes con integración con APIs de mapas

Aunque podemos también cambiar el tipo de mapas con los que trabajamos si le especificamos de dónde puede sacar los tiles

m <- leaflet(data=centros) %>% 
  addTiles() %>% 
  addMarkers(~lon, ~lat, popup = ~ as.character(ETIQUETA)) %>%
  addProviderTiles(providers$CartoDB.Positron)
m

Paquetes para trabajar con objetos geográficos

Hay otra manera de visualizar mapas sin utilizar imagenes (raster). Mediante polígonos que definen fronteras (paises, regiones, ciudades, etc.). Estos polígonos también se pueen utilizar para definir carreteras, accidentes geográficos, edificios, puntos de interés, etc. El formato más utilizado para trabajar con estos objetos espaciales es el de shapefiles. Otro muy extendido es GeoJSON o TopoJSON.

Hay muchos sitios para descargar los shapefiles de diferentes regiones administrativas, a diferentes niveles de resolución Por ejemplo:

Paquetes para trabajar con objetos geográficos

Hay muchos paquetes para trabajar con objetos geográficos en R. Los más importantes que vamos a ver son rgdal y sp. Por ejemplo, vamos a bajarnos los shapefiles de los barrios de Madrid y los visualizamos. Para cargar los shapefiles vamos a utilizar el paquete rgdal package

require(sf)
map_barrios <- read_sf("./data/Barrios Madrid/200001469.shp")

Paquetes para trabajar con objetos geográficos

Cuando cargamos un shapefile en R con sf una de las columnas de la tabla son los polígonos. Es importante también comprobar la proyección utilizada:

head(map_barrios,3)
## Simple feature collection with 3 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 438771 ymin: 4473158 xmax: 441379.7 ymax: 4475201
## Projected CRS: ED50 / UTM zone 30N
## # A tibble: 3 × 4
##   CODBDT GEOCODIGO DESBDT                                               geometry
##    <int> <chr>     <chr>                                           <POLYGON [m]>
## 1 797664 079011    011 Palacio     ((439930.8 4475097, 439943.1 4475086, 439956…
## 2 797665 079012    012 Embajadores ((440478.8 4474145, 440490 4474140, 440496.7…
## 3 797666 079013    013 Cortes      ((440981.7 4474699, 441019.1 4474673, 441023…
st_crs(map_barrios)
## Coordinate Reference System:
##   User input: ED50 / UTM zone 30N 
##   wkt:
## PROJCRS["ED50 / UTM zone 30N",
##     BASEGEOGCRS["ED50",
##         DATUM["European Datum 1950",
##             ELLIPSOID["International 1924",6378388,297,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4230]],
##     CONVERSION["UTM zone 30N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Europe - between 6°W and 0°W - Channel Islands (Jersey, Guernsey); France offshore; Gibraltar; Ireland offshore; Norway including Svalbard - offshore; Spain - onshore; United Kingdom - UKCS offshore."],
##         BBOX[35.26,-6,80.49,0.01]],
##     ID["EPSG",23030]]

Paquetes para trabajar con objetos geográficos

Hay muchas funciones para trabajar, analizar y modificar estos objetos. Por ejemplo

dim(map_barrios) # nos dice cuantas `features` hay. En este caso polígonos
## [1] 128   4
st_area(map_barrios) # nos da el área de cada uno de los polígonos
## Units: [m^2]
##   [1]   1461419.6   1026606.4    591318.7    740494.4    951002.5    445247.6
##   [7]    967645.5   1075925.9    567725.8   1392299.7   1056584.8    646806.5
##  [13]    795475.6    744061.6    620354.6   1032389.7    491627.1   1897360.6
##  [19]    642638.4    870370.9    771045.1    854034.6   1603158.7    526446.9
##  [25]    774435.9   1707104.8   1035494.0    761851.7   1714637.4   1787306.4
##  [31]   2164233.2    713526.0   1185339.9    705792.0    999170.0   1163937.9
##  [37]    604598.0    508225.5    579176.4    611027.2    937723.6    976850.0
##  [43]   1067818.3 187812697.2   1496684.2   2860411.4   1363616.2   2156862.7
##  [49]   8984663.5   6906049.7  26365060.7  17454584.6    757918.2  14180481.8
##  [55]   1406882.6   3292882.8   3558004.3   5833273.1   1292174.4   1389022.4
##  [61]   1684427.5   2839859.4   9206193.8   5421319.4   3605761.4    664786.2
##  [67]   1109593.1   1901910.0   1598661.5   1606215.9   5625327.4   1570519.9
##  [73]   1313311.9   1422619.3   1470789.0    775134.7    902190.0    768550.5
##  [79]   1101274.8   5990050.9   1072588.5   1720319.2   3071568.1   1248388.2
##  [85]   1857568.9   1011934.8    738124.3   1793685.3    996243.3    970546.8
##  [91]    591806.0   3191455.8   2325994.8    722078.1    888438.5   1052470.5
##  [97]   1015481.2    566551.4    245188.1   1427627.6   1143402.8   3146667.0
## [103]   2523248.8   2622101.6   1230948.2  16727328.1   9082299.0   1093975.5
## [109]   6364719.4   1524663.8   2112397.5  49694871.0   2062994.9  32561020.4
## [115]   2424385.0   2271836.7    548334.9    367851.2   1304313.1   9258990.6
## [121]   5004149.1   1617366.0   1880616.5   1983858.2  25095226.9    628745.8
## [127]   9568166.0   4660164.0

Por otro lado st_centroid(map_barrios) nos devuelve las coordenadas de los centros de cada polígono

cc <- st_centroid(map_barrios)
head(st_geometry(cc),2)
## Geometry set for 2 features 
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 439528.8 ymin: 4473624 xmax: 440497.1 ymax: 4474317
## Projected CRS: ED50 / UTM zone 30N

Como vemos las coordenadas no están en lat, lon. Eso es debido a que el mapa utiliza una proyección de datos geográficos diferente a la habitual

st_crs(map_barrios)
## Coordinate Reference System:
##   User input: ED50 / UTM zone 30N 
##   wkt:
## PROJCRS["ED50 / UTM zone 30N",
##     BASEGEOGCRS["ED50",
##         DATUM["European Datum 1950",
##             ELLIPSOID["International 1924",6378388,297,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4230]],
##     CONVERSION["UTM zone 30N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Europe - between 6°W and 0°W - Channel Islands (Jersey, Guernsey); France offshore; Gibraltar; Ireland offshore; Norway including Svalbard - offshore; Spain - onshore; United Kingdom - UKCS offshore."],
##         BBOX[35.26,-6,80.49,0.01]],
##     ID["EPSG",23030]]

Paquetes para trabajar con objetos geográficos

¿Y qué es una proyección? Básicamente una proyección es un sistema de coordenadas que convierte una posición sobre el globo terráqueo en dos números (longitud y latitud). El problema es que hay varias maneras de hacerlo

  • Proyección UTM
    • En vez de la tradicional de longitud y latitud da el punto en metros respecto a líneas (paralelos y meridianos) definidos sobre la superficie.
  • Proyección WGS84
    • Es un tipo de proyección esférica, donde cada punto se da en longitud y latitud que indican dónde está sobre la esfera. El GS84 es un estándar de cómo se define la esfera terráquea. se refiere a la diferente definición de esfera. En este caso es un elipsoide
    • Esta es la tradicional

1984 World Geodetic System revision

Paquetes para trabajar con objetos geográficos

Vamos a cambiar la proyección UTM que tiene nuestro shapefile a uno con WGS84

mapb <- st_transform(map_barrios,"+proj=longlat +datum=WGS84")
mapb <- st_transform(map_barrios,4326)

y ahora como vemos las coordenadas están en longitud/latitud

cc <- st_centroid(mapb)
head(st_geometry(mapb),2)
## Geometry set for 2 features 
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -3.722965 ymin: 40.40503 xmax: -3.692199 ymax: 40.42339
## Geodetic CRS:  WGS 84

Paquetes para trabajar con objetos geográficos

Podemos visualizar todos los polígonos o solo parte de ellos

plot(mapb["DESBDT"])

plot(mapb[mapb$DESBDT %in% c("011 Palacio","015 Universidad"),"DESBDT"])

Paquetes para trabajar con objetos geográficos

Podemos poner juntos imágenes raster con puntos (ya lo hemos visto) y también los polígonos. Por ejemplo, usando ggmap podemos poner los polígonos de los distritos de Madrid encima de un mapa raster. Para ello tenemos que convertir el objeto SpatialDataFrame a un data.frameutilizando el comando fortify

ggmap(mm)+geom_sf(data=mapb,inherit.aes = F,fill="grey",size=.2,alpha=.5,color="green")

Paquetes para trabajar con objetos geográficos

Lo mismo podemos hacer con leaflet

leaflet(mapb) %>% addTiles() %>% addPolygons(weight=2,color="green",fillColor="gray",fillOpacity = 0.5)

Paquetes para trabajar con objetos geográficos

Con leaflet es muy facil también mostrar diferentes capas. Por ejemplo vamos a mostrar los distritos y los centros comerciales

leaflet(mapb) %>% addTiles() %>%
  addPolygons(weight=2,color="green",fillColor="gray",fillOpacity = 0.5) %>%
  addCircles(data=centros, lng= ~lon, lat=~lat,popup = ~ ETIQUETA)

Choropleths

Una de las visualizacións de información geográfica más utilizada es dibujar los polígonos con colores proporcionales a una variable (choropleths). Por ejemplo, vamos a pintar los barrios con color proporcional a la población en cada uno.

Primero nos bajamos la población de los barrios de Madrid de
http://www-2.munimadrid.es/TSE6/control/seleccionDatosBarrio

pob_barrios <- read.csv("./data/Madrid-Poblacion_Barrios.csv",
                        colClasses=c("character","character","integer","character","character"))

Añadimos a la tabla del mapa una columna que sea la población (uniendo por el GEOCODIGO)

mapb <- merge(mapb,pob_barrios,by="GEOCODIGO")
head(mapb)
## Simple feature collection with 6 features and 7 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -3.722965 ymin: 40.40503 xmax: -3.690376 ymax: 40.43061
## Geodetic CRS:  WGS 84
##   GEOCODIGO CODBDT          DESBDT Distrito      Barrio Población     Fecha
## 1    079011 797664     011 Palacio   CENTRO     PALACIO     22471 10/1/2013
## 2    079012 797665 012 Embajadores   CENTRO EMBAJADORES     46264 10/1/2013
## 3    079013 797666      013 Cortes   CENTRO      CORTES     10717 10/1/2013
## 4    079014 797667    014 Justicia   CENTRO    JUSTICIA     16847 10/1/2013
## 5    079015 797668 015 Universidad   CENTRO UNIVERSIDAD     31714 10/1/2013
## 6    079016 797669         016 Sol   CENTRO         SOL      7577 10/1/2013
##                         geometry
## 1 POLYGON ((-3.709384 40.4224...
## 2 POLYGON ((-3.702836 40.4139...
## 3 POLYGON ((-3.696961 40.4189...
## 4 POLYGON ((-3.699416 40.4284...
## 5 POLYGON ((-3.712235 40.4302...
## 6 POLYGON ((-3.70416 40.42019...

Choropleths

Podemos mostrar el choropleth utilizando ggplot y ggmap

Choropleths

Como siempre, mucho mejor con el paquete leaflet. El paquete `leaflet’ también tiene opciones para hacer las paletas de colores de manera automática (https://rstudio.github.io/leaflet/colors.html)

pal <- colorNumeric(palett="magma",domain=mapb$Población)
m <- leaflet(data = mapb) %>% 
  addPolygons(fillColor= ~ pal(Población),weight=1,fillOpacity = 0.8,label=~as.character(DESBDT))
m

Choropleths

Y hay muchas opciones para hacer más interesantes el resultado final (https://rstudio.github.io/leaflet/choropleths.html). Por ejemplo, podemos añadir una leyenda para informar sobre los valores

pal <- colorNumeric(palett="magma",domain=log(mapb$Población))
m <- leaflet(data = mapb) %>% 
  addPolygons(fillColor= ~ pal(log(Población)),weight=1,fillOpacity = 0.8,popup=~as.character(DESBDT)) %>% addLegend(pal=pal,values=~log(Población))
m

Geolocalizar puntos

Una de las operaciones más usuales en datos espaciales es encontrar el polígono (polígonos) donde se encuentra un punto. Para eso utilizamos la función st_within en el paquete sf

Por ejemplo: ¿en qué barrio está el Bernabeu y el Calderón?

gps_Bernabeu <- c(-3.687909,40.4544)
gps_Calderon <- c(-3.720589,40.40671)

Transformamos estos puntos a objetos de sf

pto_matrix <- rbind(gps_Bernabeu,gps_Calderon)
estadios <- st_as_sf(data.frame(pto_matrix),coords=c("X1","X2"),crs=4326)

Finalmente encontramos el índice (número de polígono) en el que están los puntos y el barrio

id_polygon <- st_within(estadios,mapb)
mapb[as.numeric(id_polygon),]
## Simple feature collection with 2 features and 7 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -3.722965 ymin: 40.39966 xmax: -3.663164 ymax: 40.46004
## Geodetic CRS:  WGS 84
##    GEOCODIGO CODBDT             DESBDT   Distrito         Barrio Población
## 29    079054 797692 054 Hispanoamérica  CHAMARTIN HISPANOAMERICA     31374
## 7     079021 797670       021 Imperial ARGANZUELA       IMPERIAL     22978
##        Fecha                       geometry
## 29 10/1/2013 POLYGON ((-3.664455 40.4588...
## 7  10/1/2013 POLYGON ((-3.71465 40.40017...

Calcular distancias y áreas.

Otra de las operaciones más habituales en datos geo-espaciales es medir distancias de un punto a otro o el área de un polígono.

Para medir la distancia entre dos puntos, utilizaremos la función spDists o spDistsN2 del paquete sp. Por ejemplo, esta es la matriz de distancias entre el Bernabeu y el Calderon

st_distance(estadios) #calcula todas las distancias entre todas las filas de pto_matrix
## Units: [m]
##          [,1]     [,2]
## [1,]    0.000 5980.953
## [2,] 5980.953    0.000

o entre todas las filas y una de ellas en particular

st_distance(estadios,estadios[1,])
## Units: [m]
##          [,1]
## [1,]    0.000
## [2,] 5980.953

Calcular distancias y áreas.

Lo mismo para las áreas. Podemos utilizar la función area del paquete raster que da el área en metros cuadrados. Por ejemplo, esto es lo que ocupa en km2 el barrio de Los Jerónimos

st_area(mapb[grep("Jerónimos",mapb$DESBDT),],) /1000000
## 1.896524 [m^2]

mientras que él área del barrio del Pardo es

st_area(mapb[grep("Pardo",mapb$DESBDT),]) /1000000
## 187.7198 [m^2]

Ejercicio 1

  • Repetir el proceso de bajada de los shapefiles y de la población para las secciones censales de la comunidad de Madrid.

  • Mostrar en una visualización tipo choropleth la población de cada sección censal.

  • Mostrar en otra visualización tipo choropleth el porcentaje de población extranjera en cada sección censal.

Ejercicio 1

Para obtener los datos de las secciones censales tenemos que ir a la página del INE: http://www.ine.es/censos2011_datos/cen11_datos_resultados_seccen.htm

Donde tenemos dos tipos de archivos:

  • Indicadores para secciones censales (en formato CSV mejor). Mirad también el archivo de “Relación de Indicadores” para saber cuáles son.
  • Contorno de las secciones censales en formato shapefile.

Ejercicio 1

Cargamos los datos de las secciones censales de la comunidad 13 (Madrid)

indicadores <- read.table("./data/indicadores_seccion_censal_csv/C2011_ccaa13_Indicadores.csv",sep=",",header=T)

Y los shapefiles de las secciones censales

map_censales <- read_sf("./data/cartografia_censo2011_nacional/SECC_CPV_E_20111101_01_R_INE.shp",stringsAsFactors=F)

Ejercicio 1

Como antes las secciones censales están en otra proyección

st_crs(map_censales)
## Coordinate Reference System:
##   User input: ETRS89 / UTM zone 30N 
##   wkt:
## PROJCRS["ETRS89 / UTM zone 30N",
##     BASEGEOGCRS["ETRS89",
##         ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
##             MEMBER["European Terrestrial Reference Frame 1989"],
##             MEMBER["European Terrestrial Reference Frame 1990"],
##             MEMBER["European Terrestrial Reference Frame 1991"],
##             MEMBER["European Terrestrial Reference Frame 1992"],
##             MEMBER["European Terrestrial Reference Frame 1993"],
##             MEMBER["European Terrestrial Reference Frame 1994"],
##             MEMBER["European Terrestrial Reference Frame 1996"],
##             MEMBER["European Terrestrial Reference Frame 1997"],
##             MEMBER["European Terrestrial Reference Frame 2000"],
##             MEMBER["European Terrestrial Reference Frame 2005"],
##             MEMBER["European Terrestrial Reference Frame 2014"],
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]],
##             ENSEMBLEACCURACY[0.1]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4258]],
##     CONVERSION["UTM zone 30N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Europe between 6°W and 0°W: Faroe Islands offshore; Ireland - offshore; Jan Mayen - offshore; Norway including Svalbard - offshore; Spain - onshore and offshore."],
##         BBOX[35.26,-6,80.49,0.01]],
##     ID["EPSG",25830]]

Asi que la cambiamos a WGS84

mapc <- st_transform(map_censales,4326)

y os quedamos solo con las de la Comunidad de Madrid

mapc <- mapc[mapc$CCA == 13,]

Ejercicio 1

Para unir los datos de los indicadores y de los shapefiles vamos a utilizar el indicador único CUSEC. Pero el archivo de los indicadores no lo tiene. Se lo ponemos

require(stringr)
poblacion <- data.frame(CUSEC = paste(
                          str_pad(indicadores$cpro, 2, pad = "0"),
                          str_pad(indicadores$cmun, 3, pad = "0"),
                          str_pad(indicadores$dist, 2, pad = "0"),
                          str_pad(indicadores$secc, 3, pad = "0"), sep=""),
                          poblacion = indicadores$t1_1)

Y los juntamos

mapc <- merge(mapc,poblacion,by="CUSEC")

Ejercicio 1

Hacemos un choropleth con la población por seccion censal

pal <- colorNumeric("magma",domain=log(mapc$poblacion))
leaflet(mapc) %>% addTiles() %>%
  addPolygons(fillColor = ~pal(log(poblacion)),
              fillOpacity = 0.8,weight=0.1)